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ABSTRACT 


Effects  of  a  turbulent  gas  flow  on  the  stability  of  a 
glow  discharge  are  investigated.   While  it  is  known  that 
turbulence  affects  the  stability  of  glow  discharges,  the 
mechanism  is  not  clear.   Primarily,  the  problem  lies  in  the 
very  long  characteristic  times  of  turbulence  phenomena  when 
compared  to  the  glow  discharge  instability  times. 

A  mathematical  model  is  developed  which  is  solved  numer- 
ically for  the  ambipolar  diffusion  and  an  overall  energy 
equation  in  an  unsteady,  cylindrical  coordinate  system. 
Strong  perturbations  of  the  electric  field  are  introduced 
which  disturb  an  otherwise  stable  configuration  and  the 
effects  of  turbulence  on  the  time  evolution  of  the  perturba- 
tion are  observed. 

It  is  shown  that  the  modification  of  the  ambipolar 
diffusion  coefficient  and  the  thermal  conductivity  is  a 
reasonably  sufficient  model  to  introduce  the  turbulence 
effects.   It  is  found  that  the  charged  particle  density  is 
effectively  the  sole  source  of  heating  in  the  gaseous  environ- 
ment of  the  discharge.   It  is  then  shown  that  turbulence 
acts  to  suppress  the  temperature  instability  introduced  by 
the  discharge  streamer. 
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I.   INTRODUCTION 

If  a  gas  such  as  nitrogen  or  argon,  or  a  mixture  of 
gases  such  as  air,  is  subjected  to  an  electrical  potential, 
a  sequence  of  phenomena  take  place  as  a  result  of  a  continu- 
ously increasing  applied  voltage  and/or  current-   Initially, 
for  a  small  voltage,  a  very  weak  current  on  the  order  of  0.1 
picoamperes  develops.   This  current  will  increase  slightly 
with  increasing  voltage  until  all  the  naturally  occurring 
charge  carriers  are  being  utilized.   A  current  saturation 
condition  now  prevails  until  the  applied  voltage  reaches 
about  6  kV  in  air  at  which  time  a  further  voltage  rise 
rapidly  increases  the  current  flow  to  about  one  picoampere. 
This  is  referred  to  as  the  Townsend  Discharge  Region  and  it 
is  due  to  new  charge  carriers  being  created  by  collisions 
between  the  highly  energetic  electrons  (fast  moving)  and  the 
neutral  (slow  moving)  atoms.   This  process  can  have  several 
effects  on  the  individual  gas  molecules.   Essentially,  it 
can  excite  them,  dissociate  them,  and/or  ionize  them.   This 
ionization  effect  is  higher  for  lower  gas  pressures  since 
the  farther  the  electron  travels  between  collisions  the 
higher  will  be  its  kinetic  energy  which  is  imparted  by  the 
acceleration  produced  by  the  applied  electric  field.   When 
the  electron's  energy  is  greater  than  that  required  for 
ionization  of  the  neutral  gas  particles  the  collisions  will 
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produce  additional  electrons  which  in  turn  gain  energy, 
produce  collisions,  and  therefore,  increase  the  total  number 
of  electrons  present.   This  process  is  called  an  electron 
avalanche . 

From  this  point,  slight  increases  in  voltage  produce 
current  increases  of  several  orders  of  magnitude  at  which 
time  the  current  becomes  self  sustaining.   Reducing  the 
external  circuit  resistance  now  increases  the  current  flow 
and  a  steady  glow  is  produced  in  the  gas.   This  glow- 
discharge  is  the  familiar  phenomenon  which  occurs  in  neon 
lights  and  has  other  practical  applications  such  as  a 
pumping  scheme  for  gas  lasers. 

Finally,  when  the  current  is  increased  still  further,  a 
highly  ionized  channel  is  formed  which  collapses  the  glow. 
If  there  is  sufficient  current  carrying  capacity  in  the 
external  circuit  a  steady  current  is  established  which 
represents  a  plasma  state  called  an  arc. 

A  full  understanding  of  the  arcing  mechanism  is  diffi- 
cult due  to  the  extremely  short  formation  times  involved. 
However,  since  delaying  the  point  of  transition  to  arcing  is 
a  highly  practical  problem,  a  method  for  modeling  it  is 
needed.   A  blend  of  theory  and  experiment  is  usually  re- 
quired.  Such  a  model  could  then  be  used  to  study  the 
effects  of  the  arc  on  and  by  its  environment. 

It  has  been  observed  that  this  glow-to-arc  transition 
can  be  delayed  by  utilizing  gas  dynamic  stabilization 
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techniques  in  flowing  gas  systems  [Ref.  1].   The  purpose  of 
this  work  is  to  establish  the  extent  of  fluid  dynamic 
stabilization  that  results  from  the  introduction  of  varying 
levels  of  fluid  dynamic  turbulence  into  the  flow.   It  will 
be  shown  that  this  method  has  beneficial  effects  on  the 
stability  of  the  glow  discharge  process  and  on  the  temperature 
distribution  in  the  gas  produced  by  the  presence  of  discharge 
streamers.   For  further  discussion  and  background  information 
relating  to  ionization  and  the  breakdown  processes  refer  to 
any  standard  reference  in  the  field  or  to  Wallace,  R.J. 
[Ref.  2]. 

An  introductory  skill  level  in  atomic  physics  plus 
fundamental  proficiencies  in  mechanics,  electrical  princi- 
ples and  numerical  analysis  should  provide  the  reader  with 
adequate  preparation  for  understanding  this  material. 

The  principles  discussed  here  have  many  practical  engi- 
neering applications.   For  maximum  power  output  of  an  elec- 
trically pumped  gas  laser,  such  as  a  carbon  dioxide  laser, 
the  glow  discharge  must  be  held  essentially  at  the  limit  of 
arcing.   Since  turbulent  flows  tend  to  delay  this  glow-to-arc 
transition,  the  laser  can  be  allowed  to  operate  at  either 
the  highest  possible  power  level  or  with  a  greater  margin  of 
safety  at  lower  power  levels.   Another  application  where  a 
stabilized  high  frequency  electrical  spark  discharge  could 
be  used  as  a  means  for  injecting  heat  into  a  propulsive  duct 
(like  a  pulsejet)  has  been  suggested  by  C.E.  Tharratt  [Ref.  3], 
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This  type  of  a  device  is  theoretically  capable  of  unlimited 
thrust.   It  would  also  have  a  low  specific  fuel  consumption 
and  could  operate  outside  the  audible  range.   The  possibility 
of  nearly  silent  operation  would  be  highly  advantageous  as 
an  aircraft  propulsion  device  for  the  future. 

Two  major  cases  will  be  discussed.   The  first  case 
concerns  the  behavior  of  a  highly  transient  discharge  in  a 
non-ionized  gaseous  environment  such  as  the  lightning  dis- 
charge.  This  first  case  concerns  a  pre-glow  discharge  type 
phenomenon  since  the  streamer  propagates  in  a  non-ionized 
medium.   However,  it  serves  to  establish  the  effects  of 
turbulence  in  a  simpler  situation  in  order  to  form  a  starting 
place  for  this  research.   The  second  case  will  look  at  dis- 
charge streamers  in  a  partially  ionized  gas  which  might  be 
used  for  the  pumping  mechanism  of  a  gas  laser.   The  method 
of  approach  to  these  problems  will  be  developed  in  the  next 
chapter  where  a  mathematical  model  will  be  presented.   The 
governing  partial  differential  equations  for  the  model  pre- 
sented will  be  solved  by  numerical  integration  techniques 
utilizing  a  computer.   A  working  Fortran  program  is  presented 
which  can  be  used  to  investigate  a  wide  variety  of  the  aspects 
of  this  problem.   Significant  results  are  then  presented  in 
graphical  format  to  demonstrate  the  stabilizing  features  of 
the  turbulent  flow. 
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II.   ANALYSIS 

A.   THE  CONTINUITY  EQUATION 

Highly  ionized  channels  which  develop  at  the  final 
stages  of  a  glow  discharge  just  prior  to  arc  breakdown  are 
called  streamers.   The  transition  to  arcing  is  then  often 
due  to  streamer-breakdown.   It  should  be  noted  here  that  the 
type  of  streamers  to  be  analyzed  are  the  non-self-sustaining 
discharge  streamers  that  precede  the  arc,  rather  than  the 
electron  avalanche  (onset  streamers)  type  that  occur  earlier 
in  the  glow  discharge  phenomena. 

Streamers  are  seen  to  issue  from  either  electrode  and 
appear  always  to  precede  the  breakdown.   These  streamers 
propagate  virtually  instantaneously  for  our  purposes  as  an 
ultrafast  ionization  wave.   They  consist  of  a  partially 
ionized  plasma  column  for  which  the  conductivity  is  dependent 
directly  on  the  charged  particle  density  within  the  streamer 
[Ref.  4].   It  will  be  shown  that  the  collisions  made  by  the 
charged  particles  within  these  streamers  are  the  main  source 
of  heating  of  the  surrounding  gas.   The  value  of  this 
charged  particle  density  may  be  determined  by  the  applica- 
tion of  the  following  continuity  equation  to  the  breakdown 
process . 


9n      a  9  ,   9n>    r7  9n       2       3  ,~  ,  . 

xr  =   —  tt—  (r  T-)    -  U  t—  -  a^n   -  a-,n  (2.1) 

9t      r  9r    9r       9r     2       3 


15 


Equation  (2.1)  gives  the  time  rate  of  change  of  the  charged 
particle  density,  n,  as  a  function  of  radial  position,  r, 
and  time,  t,  for  the  decaying  streamer.   It  is  written  in 
cylindrical  coordinates  due  to  the  natural  form  of  the 
streamer.   A  schematic  drawing  of  the  streamer  geometry  is 
depicted  in  Figure  2.1.   The  right  hand  side  (RHS)  of 
Equation  (2.1)  consists  of  three  classes  of  terms.   The  first 
term  represents  the  effect  of  diffusion  due  to  the  high 
density  gradient  within  the  plasma  channel.   D   is  the  ambi- 

3. 

polar  diffusion  coefficient.   This  ambipolar  diffusion  results 
in  the  electrons  and  ions  diffusing  out  of  the  original 
channel  at  essentially  the  same  rate.   Briefly,  the  lighter 
electrons  having  higher  thermal  velocities  leave  the  plasma 
ahead  of  the  slower,  heavier  ions.   This  produces  a  local 
positive  space  charge  which  retards  the  loss  of  electrons  and 
accelerates  the  loss  of  the  ions.   The  net  effect  is  that  the 
ambipolar  diffusion  coefficient  is  approximately  twice  the 
diffusion  coefficient  of  the  ions  alone  [Ref .  4] . 

Next  consider  terms  3  and  4  of  the  RHS  of  Equation  (2.1) 
which  represent  two  and  three-body  recombination  respectively. 
Recombination  is  a  de- ionization  process  which  is  then 
responsible  here  for  the  reduction  in  conductivity  of  the 
partially  ionized  gas.   These  processes  occur  when  an  electron 
and  an  ion  collide  at  relatively  low  velocity  and  recombine 
to  form  a  neutral  atom.   In  order  to  conserve  linear  momen- 
tum a  third  body  must  be  present.   In  two-body  recombination 
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Figure  2.1.   Streamer  Geometry 
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this  third  body  is  a  photon  and,  therefore,  this  process  is 
sometimes  called  radiative  recombination.   For  three-body 
recombination  the  third  body  is  usually  an  ion  or  neutral 
atom  [Ref.  4].   Due  to  the  large  complexity  of  electron  loss 
mechanisms  accurate  values  for  the  recombination  coefficients 
are  difficult  to  obtain.   The  values  used  in  this  thesis 
will  be  discussed  in  the  next  subsection. 

The  second  term  on  the  RHS  is  the  convection  term.   This 
term  is  neglected  since  the  characteristic  time  of  convection 
(length/flow  velocity)  is  much  less  than  the  characteristic 
recombination  times  associated  with  terms  3  and  4  [Ref.  2]. 

While  the  convection  term  in  Equation  (2.1)  can  be 
neglected,  the  diffusion  term  cannot  since  the  diffusion 
coefficient  D   can  increase  by  about  three  orders  of  magni- 

a 

tude  between  laminar  and  turbulent  flow  [Ref.  2].   This  will 
cause  the  characteristic  diffusion  time  to  be  on  the  same 
order  of  magnitude  as  the  recombination  time  and  some  rela- 
tionship will  exist  between  diffusion  and  recombination  which 
will  affect  the  charged  particle  density  distribution,  n, 
in  space  and  time. 

Once  Equation  (2.1)  is  solved  for  n(r,t)  the  resulting 
distribution  profile  can  be  integrated  to  determine  a 
conductance  of  the  decaying  streamer  as  a  function  of  the 
streamer  radius  and  time.   This  can  be  accomplished  under 
varying  turbulence  level  conditions  by  changing  the  effec- 
tive value  of  the  diffusion  coefficient. 
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Equation  (2.1)  describes  the  lightning-discharge  type 
phenomenon  in  a  non-ionized  medium.   A  more  interesting 
situation  is  the  case  of  the  streamer  as  it  appears  in  an 
ionized-sustainer  discharge  for  the  pumping  of  a  gas  laser 
[Ref.  5].   For  this  case,  Equation  (2.1)  may  be  written  as, 


8n     a  9   ,   9n.  2       3  , 0  ~ , 

3t  "  T   17  (r  37)   =   Vin  "  a2n   "  a3n  {2'2) 


where  now  a  new  first  term  is  added  to  the  RHS  which  repre- 
sents a  production  of  electrons  and  ions.   The  ionization 
coefficient,  v.,  is  a  function  of  the  sustainer  E/N  and  can 
be  determined  from  Figure  2.2.   The  electric  field  referred 
to  here  is  an  existing  axial  electric  field,  that  is,  the 
"pumping  field"  that  maintains  the  glow.   Equations  (2.1) 
and  (2.2)  are  non-linear  partial  differential  equations 
which  will  be  solved  numerically  by  the  use  of  a  finite 
difference  approximation  scheme. 

B.   THE  ENERGY  EQUATION 

In  addition  to  the  continuity  equation,  a  suitable 
energy  equation  must  be  considered  in  order  to  describe  the 
temperature  distribution  introduced  into  the  gas  by  the  glow 
discharge  streamer.   For  thermal  stability  to  exist,  a 
balance  must  be  created  such  that  the  time  rate  of  change  of 
the  internal  energy  is  zero  because  the  loss  of  heat  through 
thermal  conductivity  and  the  addition  of  energy  from  Ohmic 
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Figure  2.2.   Relevant  Properties  of  Nitrogen 
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heating  must  balance  each  other.   The  following  equation 
represents  this  relationship. 


P  §f   =   E  •  j  -  V  •  q  (2.3) 


The  first  term  on  the  RHS  of  Equation  (2.3)  represents 
Joule  or  Ohmic  heating  where,  E,  is  the  electric  field 
strength  and,  j,  is  the  current  density.   The  second  term 
is  the  heat  flux  loss  due  to  thermal  conductivity  and 
radiation.   The  total  derivative, 


D£=9£+v^+v^+v^e-  (2  4 

Dt      3t     r  9r     6  38     z  3z  K    ' 


d  e 
reduces  to,  rr,  since  the  model  under  consideration  here 

o  t 

allows  for  conduction  only  in  the  "z"  direction.   It  is 
further  assumed  that  the  change  of  the  internal  energy  along 
the  length  of  the  streamer  is  small,  the  radiation  emitted 
primarily  perpendicular  to  the  streamer  centerline  has  a 
negligible  effect  and,  due  to  the  short  characteristic  times 
involved,  the  heating  occurs  at  an  essentially  constant 
volume.   With  these  assumptions  in  mind,  Equation  (2.3) 
becomes , 


pC   ||  -  KV2T   =  oeI  (2.5) 

VOL.  Z 


where,  p,  is  the  gas  density,  o,  is  the  electrical  conduc- 
tivity of  the  plasma  (recalling  that  j  =  oE  )  and  K  is  the 
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thermal  conductivity.   Equation  (2.5)  gives  the  temperature, 
T,  as  a  function  of  the  radius  and  time. 
The  conductivity,  a,  is  given  by, 


a   =   e  n  y  (2.6 

e  e 


where,  e,  is  the  elementary  electronic  charge  in  coulombs, 

n  ,  is  the  number  of  electrons  and,  u  ,  is  the  electron 
e  e 

mobility.   The  mobility  is  also  a  function  of  E/N  and  is 

defined  as  the  ratio  of  the  electron  drift  velocity,  v.,  to 

a 

the  applied  external  electric  field.   The  value  of  v.,  is 

d 

given  in  Figure  2.2  as  a  function  of  E/N. 

2 

It  will  be  shown  later  that  for  the  aE      term  from 

z 

Equation  (2.5) ,  E   is  essentially  constant  and  since  o  is 
directly  proportional  to  n  the  solution  of  Equation  (2.2) 
represents  the  single  source  of  heating  that  drives  the 
temperature  distribution  obtained  from  Equation  (2.5).   The 
consequence  of  this  is  that  a  mathematical  model  describing 
the  temporal  and  spatial  characteristics  of  the  streamer  can 
now  be  developed.   Furthermore,  the  effects  of  turbulence  as 
a  stabilizing  medium  can  then  be  studied  since  the  thermal 
conductivity,  K,  and  the  ambipolar  diffusion  coefficient, 
D  ,  are  the  only  parameters  affected  by  turbulence  in  Equations 
(2.1)  ,  (2.2)  and  (2.5)  . 


C.   SOLVING  THE  EQUATIONS 

An  effective  means  of  generalizing  a  problem  such  as 
this  is  to  non-dimensionalize  the  equations.   To  begin  then, 

22 


after  some  algebraic  manipulation,  Equation  (2.1)  can  be 
written  in  the  form, 


1   3n    ,        Da  ,3   *  3n.       -2    a3no  -3  ,,  - . 

— r  -  [ =-]— -(r  — )   =   -n   -  — n  (2.7) 


a„n   3t         2      ~    .^  a„ 

2  o       n  a„r   3r    3r  2 

o  2  o 

where  n  =  n/n   and  r  =  r/r  . 

o  o     o        o 

values  for  the  charged  particle  density  and  the  streamer 
radius.  The  following  additional  definitions  will  prove 
useful  in  this  model, 

2 
a_n  r 

6   =    -  °  °  (2.8) 

a 

a   E   -^-2.  (2.9) 

a2 

t   E   —  (2.10) 

The  parameter,  3/  is  an  inverse  measure  of  the  turbu- 
lence level  since  the  ambipolar  diffusion  coefficient  is 
smaller  in  laminar  flow  (large  6)  and  larger  in  a  turbulent 
flow  (small  6).   If  one  neglects  the  diffusion  term  in 
Equation  (2.1),  then  for  two-body  recombination,  direct 
integration  yields, 


n 

n   =   ,  ^   °   .  (2.11 

1  +  a~n  t 
2  o 


This  ion  loss  formula  [Ref.  2]  then  shows  that  the 
characteristic  recombination  time,  tr,  is  simply  l/a~n 
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Incorporating  the  above  definitions  allows  Equation  (2.7) 
to  be  written  in  the  non-dimensionalized  form, 


9n    19  .-  9n>       "2     "3  ._  ,_. 

—   ~    a    — (r  — )   =   "n   "  an  (2.12) 

9t    p  9r    9r 


which  lends  itself  to  a  numerical  solution  by  substituting 
the  appropriate  finite  difference  approximation  formula 
[Ref.  6]  for  the  differentials.   After  some  manipulation, 
Equation  (2.12)  can  be  approximated  by, 


n        =   {At[^(A+B)  -C]}  +  n.   .  (2.13) 

1  /  J  1~-L  P  1  ,  J 


where, 


A   =   -^(n.  .  .  -  2n.  .  +  n.  ,  .)  (2.14) 

A-2   i+l,:      1,3     i-l,D 


B   =   — i—  (n\  ,  .  -  n.  ,  .)  (2.15) 

~~  .  -   l+l , i  l-l , i 

2r Ar       J         J 


C   =   (n2  .  +  an3  .)  (2.16) 

i,d    i,: 


In  this  scheme,  ' i'  represents  position  perpendicular  to 
the  streamer  centerline  and  'j'  represents  relative  position 
in  time  after  some  initial  condition.   Equations  (2.13) 
through  (2.16)  form  the  basic  algorithm  used  by  the  Fortran 
program  discussed  in  Appendix  A. 
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The  accuracy  of  the  charged  particle  density  obtained  by 
this  type  of  differencing  scheme  depends  on  the  step  sizes 
At  and  Ar .   The  steps  must  be  smaller  than  the  smallest 
characteristic  time  of  the  phenomena  described.   Therefore, 
At  is  taken  as  0.1  or  in  real  time  xR/10  [Ref.  6] .   Primarily 
for  programming  convenience,  Ar  is  also  taken  as  0.1.   Addi- 
tionally, the  value  of  6  is  limited  by  the  following  stability 
criterion  [Ref.  2], 


S   >  ^  (2.17) 

Ar 


Returning  to  Equation  (2.2),  it  is  seen  that  the  addition 
of  the  ionization  term  to  the  RHS  requires  a  modification 
to  Equation  (2.16)  which  now  becomes, 


C,   =   [n2  .  +  an3  .  -  n.  .]  (2.18) 

1        i,D      1,3     1,3 


Equation  (2.5)  can  be  non-dimensionalized  in  a  similar 
manner  by  writing  it  in  the  form, 

^2 

To  9T     r  KTo  1  1  9   ~  9T,      aEz  ,0  1Q, 

T —  "  [ 2]  ~  ~[r  ~]       =      ^C~  (2*19) 

A  3t_     pC  r    r  3r    9r       p^v 
A     M  v  o 

where  t   is  defined  as, 

T  pC 
xA   =   —2-  (2.20) 

aoEz 
o 
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A  check  of  Equation  (2.20)  reveals  that  t   has  the 
units  of  seconds  and,  therefore,  is  a  characteristic  time 
for  the  energy  equation.   if  one  now  defines  the  following 
additional  relations , 


tA   E   -L  (2.21) 


T   =  ^-  (2.22) 

o 

0Ez  ro 

BA   E  ~W^  <2-231 

o 


then,  it  can  be  shown  that,  Equation  (2.19)  becomes  the 
following  non-dimensional  relation, 

4 ^^(r  li,   =   -f(-^-)  (2.24) 

3tA    6Ar  3r    dr  ao   E^ 

o 

where,  T  ,  is  the  initial  temperature,  which  for  this  work, 
is  taken  as  300°  Kelvin.   It  should  be  clear  here  that,  3  , 
is  now  a  non-dimensional  parameter  whose  value  depends  only 
on  the  thermal  conductivity,  K.   This  parameter  then  serves 
the  same  function  for  the  energy  equation  as  does  3  for  the 
continuity  equation,  that  is,  it  allows  the  influence  of 
turbulence  to  be  introduced  into  the  model.   The  value  of 
the  electric  field  strength,  E  ,  from  the  RHS  of  Equation 
(2.24)  can  be  estimated  by  assuming  that  the  initial  value 
of  the  externally  applied  field,  E  ,  is  given  by  [Ref .  8] , 
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E    =   E  +  ^p-  (2.25) 


A  reasonable  value  for  Q./1    might  be  on  the  order  of 
10,000  0,/m,    for  which  Equation  (2.25)  becomes, 


E   •=   E[l  +  ^(10,000) ]  (2.26) 

O  £j 


where  I/E  for  the  streamer  can  be  shown  to  be, 


e(y  N)    r 
|  =   f—   /   2^rne  dr  (2.27) 


and, 


N   =   j^r  (2.2: 


"N"  is  the  total  number  of  particles  present  in  the  gas  at  a 
pressure,  P,  and  a  temperature,  T.   The  Boltzman  constant  is 
k  =  1.38E-23  J/°K.   The  integral  in  Equation  (2.27)  is  a 
measure  of  the  conductance  of  the  streamer  and  is  evaluated 
in  Chapter  III  by  integrating  the  charged  particle  density 

distribution  obtained  from  Equation  (2.1)  .   The  initial 

2 
value  of  this  integral  is  simply  Tin  r  .   Using  estimates 

of  n   and  r   to  be  3.]25E18  m    and  4.0E-4  m  respectively 
o      o 

[Ref.  7],  Equation  (2.26)  can  be  solved  for  E   becoming, 
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EZQ      ■       l+0.g006       =       Eo  <2'29 


where  the  value  of  y  N  in  Equation  (2.27)  can  be  obtained 

e       ^ 

from  reference  tables  to  be  5.257E24  (v-m-s)~   for  nitrogen. 
Results  to  be  discussed  in  Chapter  III  show  that  the  value 
of  the  integral  drops  rapidly  to  a  small  and  nearly  steady 
value  (on  the  order  of  0.01)  such  that  the  denominator  of 
Equation  (2.29)  approaches  unity  in  very  short  times. 

The  results  of  the  previous  paragraph  then  show  that  the 

RHS  of  Equation  (2.24)  is  simply  n  since  o/o      is  (en  jj  /en  y  ) 

2   2 

and  E  /E    is  approximately  one.   So  now  finally,  Equation 

o 
(2.24)  can  be  written  as, 


^ —  —  [r  ^]   =   n  (2.30) 

3tA    6Ar  3r    9r 


Equation  (2.30)  then  shows  that  the  temperature  distri- 
bution due  to  the  presence  of  the  streamer  is  affected  only 
by  the  charged  particle  density. 

As  for  the  continuity  equations,  Equation  (2.30)  can  be 
solved  by  introducing  the  appropriate  finite  difference 
formula  as  follows  (after  some  rearranging  of  terms) , 


?i,j  +  l   ■   AT[A(A+B)  +  n.;j]  +  i.r.  (2.31 


where, 
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A   =   -^r[T\  ,   .  -  2T.   .  +  T.  ,   .]  (2.32 

A^2   l+l,:      1,:     1-1,: 


and, 


B   =   — —  [T.  ,  .  -  T.  .  .]  (2.33) 

2rAr   1+1'^     1"1'3 


The  Fortran  program  discussed  in  Appendix  A  solves 
Equation  (2.31)  simultaneously  with  Equation  (2.1)  or  (2.2) 
depending  on  the  case  being  considered  by  using  the  new 
value  of  n  calculated  from  each  step  as  an  input  to  the 
calculation  of  the  new  temperature  value.   The  time  step  for 
the  temperature  calculation  is  adjusted  such  that  the  real 
time  step  for  the  two  equations  is  equal. 

D.   MODELING  THE  STREAMER 

The  results  of  Chapter  III  show  that  appropriate  selec- 
tion of  the  controlling  parameters  for  Equations  (2.2)  and 
(2.30)  will  produce  a  stable  system  in  space  and  time.   Then 
the  question  is :   what  effect  does  the  momentary  appearance 
of  a  streamer  have  on  this  stable  situation?   If  a  local 
perturbation  of  the  electric  field  occurs,  then  the  ioniza- 
tion coefficient  will  momentarily  be  raised  which  will  in 
turn  raise  the  charged  particle  density  in  the  region  of  the 
perturbation.   Figure  2.2  provides  the  data  for  estimating 
the  magnitude  of  this  perturbation  by  noting  the  relative 
increase  in  E/N  for  a  corresponding  increase  in  v./N.   The 
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ratio  of  these  values  can  then  be  used  to  determine  a 
simple  multiplier  (labelled  Z  and  Z  )  for  the  n  terms 
appearing  on  the  RHS  of  Equations  (2.2)  and  (2.30)  to  simu- 
late the  effect  of  the  streamer.   Therefore,  one  can  write, 


3n          13/1  3nN                 ~         ~2           ~3                                         ,„    ... 

— r  — <r  — )       =      Zn    -    n      -    an                                             (2.34) 

9t        3Ar  3r  3r 
and 


31 L.    U  IS,       =       z    ;  (2.35 

3tA        Sftr  3r 

where  Z  and  Z   are  the  multipliers  discussed.   For  the 
finite  difference  approximations,  these  two  additions  are 
easily  handled  for  Equations  (2.2)  and  (2.30)  respectively, 
as  follows. 


Cn      =       [n   .  -  anJ  .  -  Z.n.  . ]  (2.36) 

2        1,3      i,:     l  1,3 


T.   .  ,   =   At[^-(A+B)  +  Z  n.   .]  +  T.   .  (2.37) 

1/D  +  1         3A  A  1,3      1,3 


The  Fortran  program  incorporates  this  concept  by  utilizing 
a  separate  subroutine  called  "STREMR"  which  is  identical  to 
the  subroutine  "DENSTY"  except  for  the  addition  of  Equations 
(2.36)  and  (2.37) .   In  this  scheme  then,  the  streamer  is 
modeled  by  calling  up  this  subroutine  when  desired  and  then 
observing  the  reaction  of  the  system  over  time. 
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III.   RESULTS  AND  DISCUSSION 

A.   NON- IONIZED  ENVIRONMENT  CASE 

In  this  section  the  effects  of  turbulence  on  the 

streamer  in  a  non- ionized  environment  are  considered. 

Equation  (2.12)  applies  to  this  situation  for  which  Figure 

3.1  shows  a  typical  plot  of  the  charged  particle  density 

versus  radial  position  perpendicular  to  the  centerline  in  a 

fully  laminar  flow.   Note  that  the  abscissa  and  the  ordinate 

are  expressed  in  the  non-dimensional  coordinates  of  r  =  r/r 
^  o 

and  n  =  n/n  ,  respectively.   There  are  five  curves  plotted 

which  show  the  charged  particle  density  distribution  profile 

from  the  initial  conditions  at  time  zero  in  steps  as  specified 

in  the  legend  to  the  profile  at  3000  recombination  times  (TR) . 

"TR"  is  the  computer's  version  of  "t  "  and  will  be  used  here 

R 

in  the  text.   Now,  for  this  plot,  and  except  for  Figures  3.11 

and  3.12,  for  all  the  figures  in  this  section,  n   =  3.125E18 

-3  3 

m    and  a?  =  2.0E-12  m  /sec  [Ref.  7].   These  values  then 

equate  to  a  recombination  time  of  0.16  microsecond,  and  so, 
1000TR  is  equivalent  to  0.16  millisecond  in  real  time. 

In  the  laminar  case  depicted  in  Figure  3.1  it  is  clearly 
seen  that  the  centerline  charged  particle  density  drops  off 
in  time  in  accordance  with  Equation  (2.11)  which  in  non- 
dimensional  form  would  be, 


1 

n^   =   — ^—  (3.1) 

<L     1+t 
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For  example,  at  1000TR,  n-,  =  1.0E-3.   This  decay  is  due 
primarily  to  recombination  only  since,  in  this  case, 
beta  =  50,000,  and  therefore,  the  diffusion  term  becomes  very 
small.   For  longer  times,  diffusion  comes  into  play  and  some 
spreading  of  the  charge  profile  is  seen. 

Figure  3.2  shows  the  effect  of  introducing  increasing 
amounts  of  turbulence  into  the  flow.   This  plot  compares  the 
different  charged  particle  profiles  as  they  would  appear  at 
0.32  msec  for  beta  values  of  50  (fully  turbulent),  500,  5000 
and  50,000  (fully  laminar).   Note  that  the  2000TR  curve  from 
Figure  3.1  is  identical  to  the  beta  =  50,000  curve  in  Figure 
3.2.   The  effect  of  the  turbulence  is  clearly  seen  to  spread 
out  the  profile  and  to  lower  n  on  the  centerline  as  well 
[Ref .  2]  . 

As  was  mentioned  previously,  the  conductance  of  the 
decaying  streamer  can  be  determined  by  integrating  the 
charged  particle  density  over  the  radius.   This  can  be 
derived  by  "inverting"  Ohm's  law  and  introducing  the  conduc- 
tivity expression  for  a  partially  ionized  plasma  for  which 
the  current  can  now  be  expressed  as  [Ref.  4], 

2   _ 
e  n  A 

I   =   GV   =   (^f)V   =   (-— ®-)V  (3.2) 

e  e 

where  G  is  the  conductance,  or  the  reciprocal  of  the  resis- 
tance, in  the  familiar  relation,  V  =  IR.   'a'  is  the 
conductivity,  'A'  the  cross-sectional  area  and  '  2, '  is  the 
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length  of  the  streamer  (or  the  distance  between  the  elec- 
trodes in  Figure  2.1  for  this  case).   Additionally,  m   and 
v   are  the  electron  mass  and  the  electron  collision  frequency, 
respectively. 

For  the  decaying  streamer,  the  radius,  and  therefore  the 
area,  is  changing  but  the  length  is  constant.  Since  all  the 
other  factors  within  the  brackets  of  Equation  (3.2)  are  also 
constant,  'G'  can  then  be  expressed  as, 


G   =  -— (   n2irdr  (3.3) 

£m  v    ' 
e  e  0 


where  gamma  is  now  defined  as  the  value  of  the  integral  in 
the  above  equation.   Gamma  has  the  units  of  reciprocal 
meters,  and  is  determined  in  the  program,  "CHARGE"  (Appendix 
A),  by  the  subroutine,  "CNDUCT."   Gamma  is  then  directly 
proportional  to  the  conductance. 

Figure  3.3  shows  gamma  plotted  over  time  out  to  100TR. 
Note  that  gamma  drops  in  value  by  90%  in  only  about  10TR 
(1.6  microsec) .   After  this  time  the  value  levels  off  to 
about  0.04  at  100TR  and  then  beyond  the  range  of  Figure  3.3 
the  value  drops  very  slowly  to  about  0.001  at  10,000TR  (1.6 
msec) . 

Figure  3.3  is  for  fully  laminar  flow.   Introducing 
turbulence  into  the  flow  causes  a  very  slight  increase  in 
gamma  over  all  time  for  the  complete  streamer.   Figure  3.4 
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shows  this  effect  at  1000TR.   Figure  3.5  shows  the  value  of 
gamma  within  the  original  streamer  tube  channel  only.   Here 
it  is  seen  that  for  small  beta  (high  turbulence)  the  conduc- 
tance is  lower  than  in  the  laminar  flow  condition. 

Referring  back  to  Figure  3.2,  and  with  these  conductance 
characteristics  in  mind,  it  seems  rather  apparent  that  one 
'effect  of  the  turbulence  is  to  delay  the  possibility  of  a 
transition  to  arcing  since  the  conductance  is  lower  along 
the  original  channel.   Another  consideration  is  what  happens 
to  a  follow-on  streamer  appearing  after  the  charged  particle 
density  profile  of  the  first  streamer  has  died  away  for  say, 
2000TR  as  in  Figure  3.2.   In  a  laminar  flow  there  would  be 
a  distinct  path  of  least  resistance  along  the  original 
channel.   However,  in  the  turbulent  flow,  the  profile  is 
spread  out  over  a  much  larger  radius  and  then,  this  "follow- 
on  streamer"  essentially  has  no  favored  "route"  for  arcing. 

Equation  (2.30)  showed  that  the  charged  particle  density 
was  the  sole  source  of  heating  in  this  model.   Figure  3.6 

shows  the  result  of  solving  Equation  (2.30)  in  a  fully 

2 

laminar  flow  for  an  E/N  of  1.635E-17  V-cm  ,  at  one  atmosphere, 

and  300°K  for  an  electric  field  strength  of  4 . 0E4  V/m.   Note 
that  the  temperature  rise  amounts  to  less  than  about  0.0  5% 
in  1000TR  and  then  rises  very  slightly  over  time  eventually 
peaking  out  at  about  4000TR  and  then  falling  back  slowly. 
This  would  seem  to  be  a  logical  behavior  for  this  case  since 
the  charged  particle  density  that  is  heating  the  gas  is  also 
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dying  away.   The  question  is  then;  what  happens  if  one 

increases  the  electric  field? 

Figures  3.7  through  3.10  show  the  effects  of  increasing 

E   from  1.0E6  V/m  (a  typical  value  of  interest  for  lasers) 

to  8.0E6  V/m.   Note  that  $   changes  for  each  of  these  cases. 

This  is  due  to  the  fact  that  the  electron  drift  velocity  is 

a  function  of  E/N,  which  is  plotted  in  Figure  2.2.   Since  the 

drift  velocity  changes,  the  electron  mobility  changes  which, 

in  turn,  produces  a  new  conductivity  value  for  each  different 

electric  field  strength  case.   Recalling  Equations  (2.20) 

and  (2.23)  ,  one  can  see  that  each  case  has  a  different 

characteristic  time  and  a  different  B, •   For  increasing  E  , 

A  o 

xA  gets  smaller  while  S,  gets  larger. 

The  main  difficulty  with  this  is  that  the  time  step  for 
the  energy  equation  solution  must  be  adjusted  to  ensure  that 
the  continuity  equation  is  "feeding"  the  charged  particle 
density  values  into  the  energy  equation  at  the  correct  rate. 
For  the  continuity  equation,  the  non-dimensionalized  time 
step  is  0.1  (t  /10  in  real  time) .   Therefore,  the  true  time 

R 

for  one  time  step  (one  "DT"  in  the  program)  is  simply 
DT  xx_.   Similarly,  the  true  time  for  the  energy  equation  is 
DTA  XT     it  follows  then  that  the  value  needed  for  DTA  is 
simply  (DTxTr)/t  . 

Returning  to  Figures  3.7  through  3.10,  it  can  be  clearly 
seen  that  increasing  E   increases  the  temperature  rise,  as 
would  be  expected.   The  effects  of  longer  times  is  the  same 
as  the  case  for  Figure  3.6  but  requires  much  longer  computing 
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times  to  reach  the  peak  and  to  begin  falling  back.   The 

effect  of  turbulence  on  the  cases  of  interest,  say  4.0E4 

V/m  to  1.0E6  V/m,  is  to  essentially  eliminate  the  temperature 

rise . 

One  remaining  question  for  this  non-ionized  case  is  what 

effect  would  increasing  n   have?   Figure  3.11  shows  that  at 

J   o  ^ 

E   =  1.0E6  V/m  and  n   increased  by  a  factor  of  10  to  3.125E19/m 

o  '  o  2  ' 

produces  a  more  squared  off  profile  than  that  depicted  in 

Figure  3.7  suggesting  that  the  temperature  dies  rapidly 

outside  of  the  immediate  influence  of  the  streamer.   Figure 

3.12  shows  the  effect  of  another  10-fold  increase  in  N. 

Here  the  resulting  temperature  rise  drops  off  at  an  even 

steeper  rate  at  2000TR  than  for  Figures  3.7  and  3.11. 

Recalling  that  t   =  l/a~n  ,  it  is  clear  that  increasing  n 

R        Z    O  O 

decreases  the  characteristic  recombination  time  significantly 
which,  in  turn,  lowers  n  such  that  there  is  a  noticeable 
reduction  in  the  temperature  profile  outside  of  this  immedi- 
ate streamer  influence. 

B.   IONIZED  ENVIRONMENT  CASE 

In  this  section  the  effects  of  a  glow  discharge  in  an 
environment  where  a  steady  glow  exists  is  considered.   For 
this  situation,  Equation  (3.4)  is  the  applicable  non-dimensional 
form  of  the  governing  continuity  equation. 


dn  19,*  9n>      "         "2     "3  ,-  ,. 

— —  -   — -  — -(r  — )   =   n  -  n   -  an  (3.4) 

3t    Br  9r    3r 
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The  temperature  effects  are  then  determined  by  Equation 
(2.30)  where  n  is  obtained  by  solving  Equation  (3.4).   The 
first  question  in  this  case  might  be;  what  effect  does 
turbulence  have  on  the  temperature  profile  due  to  this 
steady  glow  discharge?   At  sub-atmospheric  pressure,  it  is 
known  that  the  charged  particle  density  profile  takes  on  a 
Bessel  function  distribution  [Ref.  9],   Therefore,  to  inves- 
tigate this  situation,  Equation  (2.30)  can  be  solved  inde- 
pendently from  Equation  (3.4)  by  calling  up  the  subroutine 
"DENSTY"  in  the  program  and  sending  it  a  constant,  for  all 
time,  value  for  n.   The  Bessel  function  serves  as  the  initial 
condition  profile  for  the  value  of  n.   This  distribution, 
labeled  "Nl"  in  the  program,  allows  n  to  go  to  zero  at  r  =  4.8 
Figure  3.13  shows  the  results  of  such  a  case  for  strictly 

laminar  flow.   The  values  of  the  controlling  parameters  for 

-3 
Figures  3.13  through  3.16  are:   n   =  3.125E18  m   , 

3 
a~  =  2.0E-12  m  /sec  and  E   =  1.0E6  V/m.   It  is  clearly  evident 

that  the  temperature  for  such  a  laminar  case  would  rise  to 
approximately  5400°K  (recall  that  T/T   =  1  is  equivalent  to 
300°K)  in  only  about  1000TR.   Since  molecular  nitrogen  (N  ) 
begins  to  dissociate  at  about  4000°K,  the  situation  depicted 
in  Figure  3.13  is  undesirable.   However,  introducing  turbu- 
lence into  the  flow  tends  to  significantly  reduce  this 
temperature  distribution.   Figure  3.14  shows  that  a  moderate 
degree  of  turbulence  can  be  selected  which,  essentially, 
holds  the  temperature  profile  near  the  laminar  profile  without 
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going  unstable  as  in  Figure  3.13.   Figures  3.15  and  3.16 
show  the  effects  of  increasing  the  turbulence  level  still 
further.   Note  that  in  Figure  3.15,  the  temperature  profile 
is  nearly  "frozen  out"  at  about  2400°K,  while  Figure  3.16 
is  holding  steady  at  around  1200°K. 

This  asymptotic  behavior  of  the  temperature  rise  in 
reaching  this  stable  profile  appears  to  be  related  to  the 
characteristic  time  associated  with  Equation  (2.30).   Observe 
that  for  all  three  figures  (3.14,  3.15  and  3.16)  the  tempera- 
ture rise  stabilizes  between  about  2000TR  and  3000TR.   This 
effect  will  continue  to  be  evident  in  the  remaining  cases  to 
be  discussed.   The  parameters  for  the  case  presently  being 
considered  are  such  that  the  true  time  of  2500TR  is  0.4  msec. 
This  is  approximately  4.5  times  the  characteristic  time  step 
in  seconds  for  the  energy  equation,  which  is  consistent  with 
the  concept  of  a  rise  time  in  an  asymptotic  system. 

The  final  case  to  be  discussed  here  is  perhaps  the  most 
interesting.   This  will  be  the  case  where  there  is  a  steady 
level  of  ionization  in  the  gas  and  a  perturbation  of  the 
electric  field  is  introduced  to  model  the  effects  of  the 
streamer.   For  the  steady  state  situation  at  the  centerline, 
T-r   =  0 ,  so  one  can  write  Equation  (3.5)  . 

ad,   on.       _  2       3        ._._. 

-^3?(r97\   =   °   =   Vi\  ~  a2\   "  a3n<L        (3'5) 

Neglecting  the  three-body  recombination  term  (a   is  on  the 

order  of  1.0E-34  m  /sec)  ,  and  noting  that  tt—  on  the  centerline 

'  r       dr 
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is  also  zero,  the  ionization  term  is  then  equal  to  a0n  . 
Since  Figure  2.2  gives  v . /N ,  the  charged  particle  densit 
for  the  steady  state  condition  can  be  obtained  by, 


v  .       v. 

n„  z      —     =       (~rr)  (— )  (3-6) 

<£     a2       N   a2 


where  N  is  given  by  Equation  (2.28).   Figure  3.17  shows  this 

steady  state  charged  particle  distribution  where,  for  this 

model,  the  glow  discharge  is  contained  in  a  region  of  a 

radius  equal   to  10  r  .   At  r  ,  the  boundary  condition  is  set 
^  o       o  J 

such  that  n  is  zero  at  the  "wall."   The  left  hand  boundary 

condition  assumes  that  the  slope  of  the  charged  particle 

density  is  zero.   See  Appendix  A  for  details  concerning  how 

the  program  meets  these  boundary  conditions.   Note  that  in 

Figure  3.17,  the  steady  state  condition  is  reached  in  less 

than  100TR,  and  is  constant  at  the  longer  times.   The  first 

vertical  line,  at  r/r   =  1.0,  is  the  initial  condition 

o 

distribution  which  assumes  n  =  1.0  out  to  r  =  1.0  and  zero 
beyond  this  point. 

The  steady  state  solution  for  Equation  (2.30)  is  given 
in  Equation  (3.7). 


T   =   -(B  „,  £j-+  SA  V  +  To  <3-7) 


The  boundary  conditions  are  T  =  T   at  r  =  b  and  —  =  0  at 

3r 

r  =  0.   Figure  3.18  shows  this  result  for  E   =  1.47E5  V/m. 

o 

For  this  case  (3,  =  0.2219,  and  so  at  the  centerline  r  =  0 
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and  substituting  in  the  appropriate  values  into  Equation 
(3.7)  one  obtains  T  =  6.548  which  is  confirmed  at  r  =  0  in 
Figure  3.18.   Note  that  this  steady  state  condition  is 
reached  at  approximately  2500TR  which  is  about  the  same  time 
period  required  for  stabilization  from  the  previous  case 
for  the  Bessel  function  initial  condition  profile  (Figures 
3.14,  3.15  and  3.16)  and  for  the  case  showing  the  temperature 
profile  of  a  lone  streamer  in  a  non-ionized  environment 
(Figures  3.6  through  3.12)  . 

Figures  3.19  through  3.28  show  the  results  of  introducing 
the  electric  field  perturbation  into  the  model.   Each  plot 
shows  the  situation  of  the  streamer  appearing  in  a  steady 
state  turbulent  environment  and  lasting  for  increasing 
periods  of  time  beginning  at  5  recombination  times  and  being 
incremented  by  5TR  up  to  a  final  value  of  50TR.   In  real  time, 
this  would  be  equivalent  to  0.6  msec  to  6.0  msec.   For  these 

plots,  E   =  1.47E5  V/m  which  produces  an  E/N  of  6.0E-17 

2 
V-cm  .   This  case  was  chosen  for  example  purposes  since  the 

steady  state  temperature  profile  is  convenient  at  a  centerline 
value  of  T/T   =  6.5  and  the  rise  due  to  the  perturbation 
reaches  a  maximum  of  about  13.75  (4l25°K)  at  50TR.   Each  figure 
shows  five  curves,  3  are  solid  and  2  are  dashed.   The  solid 
lines  represent  the  initial  temperature  distribution  (T  =  1.0) , 
the  steady  state  profile,  and  then,  the  profile  of  the  streamer 
at  the  end  of  its  stay  time.   The  dashed  lines  are  the  decay- 
ing, or  transient,  profiles  produced  after  the  appearance  of 
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the  streamer.   Note  that  the  centerline  temperature  drops  to 
about  half  of  the  rise  value  in  only  20TR  and  has  nearly 
reached  steady  state  in  another  100TR. 

Figures  3.29  and  3.30  show  a  case  for  a  higher  E/N  value 

2 
of  7.0E-17  V-cm  .   Here  the  steady  state  temperature  profile 

has  a  centerline  value  of  about  10.0  and  then  a  streamer  "on 

time"  of  only  8TR  sends  the  profile  above  4000°K.   Cases  for 

smaller  E/N  values  produce  lower  steady  state  temperature 

profiles.   The  streamer  effects  are  similar  except  that 

much  longer  "on  times"  are  needed  to  produce  much  of  a 

temperature  rise.   This  appears  consistent  in  the  model 

since  at  the  lower  temperature  it  would  take  a  much  larger 

perturbation  to  upset  the  steady  state  condition,  while  at 

the  higher  energy  levels  a  comparatively  smaller  perturbation 

is  all  that  is  required. 
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IV.   CONCLUSIONS  AND  RECOMMENDATIONS 

This  thesis  presents  a  mathematical  model  for  describing 
the  effects  of  varying  levels  of  turbulence  on  glow  dis- 
charge phenomena.   The  development  is  concerned  only  with 
the  aftereffects  of  the  streamer,  rather  than  its  propaga- 
tion.  The  program  that  is  presented  in  Appendix  A  is 
intended  to  allow  the  investigation  of  a  large  variety  of 
cases  of  interest  relating  to  the  glow  discharge.   The  cases 
discussed  in  the  results  chapter  are  those  that  appeared  to 
develop  a  logical  sequence  of  thought  in  attempting  to 
answer  the  question  of  how  turbulent  gas  flows  can  and  do 
affect  the  stability  of  a  glow  discharge  and  the  problem  of 
delaying  the  glow-to-arc  transition. 

The  results  clearly  show  that  the  modelled  turbulent  flow 
does  indeed  have  strong  beneficial  effects  on  the  stability 
of  the  system.   Additionally,  the  mathematical  model  success- 
fully predicts  within  a  reasonable  level  of  uncertainty  a 
region  of  stability  that  is  consistent  with  the  normal  oper- 
ating ranges  of,  say,  an  electrically  pumped  gas  laser 
[Ref .  10] .   The  model  suggests  that  the  region  for  which  the 
system  maintains  thermal  stability  is  critically  dependent 
on  E/N.   The  research  showed  that  a  steady  state  temperature 

increase  from  1200°K  to  3000°K  results  from  an  increase  in 

2  2 

E/N  of  only  5.0E-17  V-cm   to  7.0E-17  V-cm  . 
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One  important  question  still  remains  however.   That  is; 
is  there  a  particular  level  of  turbulence  which  would  have 
the  most  beneficial  effect  on  the  streamer-initiated  break- 
down problem?   This  work  suggests  that  the  best  situation  is 
simply  that  the  most  desirable  level  is  the  highest  level 
that  can  be  maintained.   However',  the  non-linear  nature  of 
the  equations  suggests  that  some  optimal  level  may  exist 
[Ref.  2].   Experimental  work  is  probably  the  best  means  of 
answering  this  question. 

The  results  presented  here  by  no  means  exhaust  the  aspects 
of  this  problem  that  could  be  investigated.   The  intent  of 
the  program  is  that  it  can  be  used  as  a  tool  to  analyze  any 
desired  case.   Some  logical  next  steps  might  be  to  expand  on 
the  streamer  modeling  aspects  to  investigate  what  happens  when 
a  follow-on  streamer  appears  before  the  initial  one  has  fully 
decayed  back  to  the  steady  state  condition.   The  "Z"  and  "Z  " 
arrays  could  be  set  up  to  have  two  or  more  streamers  appear 
at  any  desired  location.   Also,  the  values  of  "Z"  and  "Z" 
used  -in  this  model  are  somewhat  arbitrary  and  probably  have 
to  be  specified  for  each  individual  E/N  case  according  to 
experimental  observation.   More  research  into  the  actual 
magnitude  of  the  electric  field  perturbation  may  allow  the 
model  to  more  accurately  predict  the  areas  of  instability 
with  streamer  on  times  that  are  more  in  agreement  with 
experiment  [Ref.  5]. 

Mathematically,  the  analysis  of  this  problem  consists  of 
solving  the  continuity  equation  for  three  cases.   The  first 
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is  for  the  RHS  of  Equation  (3.4)  equal  to  less  than  zero. 
This  is  the  non-ionized  environment  case.   The  second  is  for 
the  RHS  equal  to  zero.   This  is  the  steady  state  case  which 
is  what  was  applied  to  the  ionized  environment  cases  discussed 
The  third  case,  which  was  not  applied  in  this  thesis,  is  for 
the  RHS  of  the  continuity  equation  to  be  greater  than  zero. 
Investigation  into  this  third  case  could  possibly  clarify 
some  aspects  of  the  problem.   And  finally,  the  energy  equation 
may  have  eigenvalue  solutions  which  might  produce  results  of 
interest  such  as  additional  regions  of  stability. 
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APPENDIX  A 
THE  CHARGE  PROGRAM 

The  Fortran  program  "CHARGE"  is  intended  to  be  a  working 
tool  for  analyzing  the  glow  discharge  problem.   It  evolved 
from  a  previously  written  program  which  solved  the  non-linear 
continuity  equation  (2.1)  [Ref .  2] .   It  is  now  arranged  in  a 
somewhat  different  form.   The  program  consists  of  a  main 
program  section  which  calls  up  any  needed  combination  of  the 
seven  subroutines  which  perform  the  desired  tasks.   As  cur- 
rently written,  it  can  compare  four  separate  cases  plus  it 
will  display  the  initial  conditions  of  either  the  charged 
particle  distribution  or  the  temperature  profile.   Individual 
test  cases  are  set  up  by  commenting  out  (placing  a  "C"  in 
column  1)  those  commands  not  desired.   After  the  main  program, 
there  are  two  sections.   The  first  provides  commands  to  print 
out  the  data  and  the  second  provides  commands  to  plot  the 
results  by  utilizing  the  "DISSPLA"  system  developed  by  the 
Integrated  Software  Systems  Corporation,  San  Diego,  California 

The  program  as  it  appears  following  this  discussion  is 
set  up  to  produce  the  plot  displayed  in  Figure  3.28.   The 
main  program  begins  by  defining  all  arrays.   For  this  exam- 
ple the  arrays  are  set  to  101  in  the  "i"  direction  which 
specifies  the  size  of  the  array  with  respect  to  radial  posi- 
tion.  Since  one  DR,  one  radial  step,  is  equal  to  r  /10, 
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position  "101"  specifies  a  distance  of  10  r   perpendicular 
from  the  centerline  of  the  streamer.   The  "j"  size  is  always 
kept  at  2  due  to  the  computing  scheme  used  in  the  subroutines 
After  setting  the  array  size,  set  the  variable  named 
"RADIUS"  equal  to  the  same  value  as  the  "i"  coordinate  above. 
In  this  case,  RADIUS  =  101.   This  determines  how  far  in  the 
radial  direction  that  the  differentiation  is  to  procede. 
Next  one  must  assign  the  desired  values  for  BETA  and  BETAA 
depending  on  the  case  at  hand.   The  BMIN ,  BMAX ,  BAMIN  and 
BAMAX  parameters  apply  when  it  is  desired  to  transition  the 
calculation  from  a  laminar  to  a  turbulent  environment.   This 
feature  is  designed  to  attempt  to  ensure  that  the  calculation 
accounts  for  the  fact  that  the  streamer  initially  appears  so 
quickly  that  the  flow  looks  laminar  to  the  streamer  no  matter 
what  the  level  of  turbulence  is.   Then,  after  a  period  of 
time,  the  charged  particle  distribution  left  behind  by  the 
streamer  has  existed  for  sufficient  time  to  now  "see"  the 
turbulent  effects.   The  time  required  for  this  transition  is 
estimated  by  noting  that  the  characteristic  diffusion  time 
for  laminar  flow  is  approximately  0.1  msec  [Ref.  2].   Depend- 
ing on  the  case  at  hand,  one  DT  is  on  the  order  of  a  few 
microseconds.   Then,  100  DT  (10TR)  is  on  the  order  of  a  few 
tenths  of  a  millisecond.   Therefore,  for  most  cases,  a 
transition  of  100  DT  is  used. 

Returning  to  the  example,  one  procedes  by  selecting  the 
desired  initial  conditions.   In  this  example,  all  of  the 
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"Nl"  values  are  commented  out  since  they  apply  to  the  Bessel 
function  initial  condition  profile.   Four  series  of  subrou- 
tine calls  are  provided  then  to  allow  running  of  four 
separate  cases.   Each  case  assigns  its  respective  results  to 
the  "VALX"  (X  =  0,1,2,3,4...)  array  for  plotting  later. 

The  initial  conditions  are  assigned  to  VALO  or  VALTO. 
Then  the  first  subroutine  call  will  see  the  initial  distribu- 
tion values  of  n  while  each  successive  call  will  see  the 
results  of  the  previous  call.   If  it  is  desired  to  recover 
the  initial  conditions,  call  either  "NZERO"  or  "TZERO"  as 
desired  and  the  current  values  will  then  be  reset  to  time 
zero . 

For  this  example  now,  "VAL1"  is  the  result  of  two  calls. 
The  first  is  for  "TRNSTY"  to  transition  from  laminar  to 
turbulent  conditions  with  time  set  to  101,  then  "DENSTY"  is 
called  and  runs  for  23,900  DT.   The  result  is  the  steady 
state  temperature  distribution  in  turbulent  flow  at  2500TR. 
VAL2  is  then  the  result  of  calling  up  the  subroutine  "STREMR" 
and  running  it  for  50  0  DT  with  BETA  and  BETAA  set  to  the 
laminar  values.   VAL3  then  results  from  calling  up  "TRNSTY" 
for  100  DT  and  "DENSTY"  for  an  additional  100  DT  at  the 
turbulent  values.   Finally,  VAL4  simply  runs  "DENSTY"  for 
another  1000  DT  in  turbulent  conditions  to  check  the  decay- 
ing temperature  profile.   Of  course,  any  combination  of  these 
subroutine  calls  and  time  values  can  be  selected  to  view  the 
situation  as  desired. 
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The  results  can  then  be  printed  out  and/or  plotted.   To 
become  familiar  with  the  "DISSPLA"  plot  system,  refer  to 
DISSPLA  Pocket  Guide  (current  with  version  9.0)  by  ISSCO 
Graphics,  10505  Sorrento  Valley  Rd ,  San  Diego,  California. 

Some  pointers  of  potential  errors  should  be  clarified. 
When  selecting  the  desired  heading  lines,  ensure  that  the 
last  entry  in  the  call  to  "HEADIN"  correctly  indicates  the 
total  number  of  heading  lines.   If  no  heading  appears,  this 
is  probably  the  reason.   Secondly,  ensure  that  the  third 
entry  in  the  call  to  "CURVE"  is  equal  to  the  "i"  value  on 
the  arrays.   This  is  the  number  of  points  (101  here)  to  be 
plotted.   And  finally,  the  second  entry  in  the  call  to 
"LEGEND"  is  the  number  of  lines  of  legend  needed,  5,  in  this 
example . 

The  subroutines  "DENSTY,"  "STREMR"  and  "TRNSTY"  all 
contain  essentially  the  same  algorithm  with  modifications  to 
perform  their  particular  function.   The  basic  idea  here  is 
that  the  calculation  procedes  from  the  first  row  of  data 
and  fills  the  second  row  with  updated  values.   Then  the  first 
row  of  old  data  is  replaced  by  the  second  row  of  new  data  by 
"SWITCH"  and  the  loop  increments  one  DT  and  procedes.   With 
this  scheme  the  program  can  be  run  as  long  as  desired  without 
using  large  values  of  memory  storage. 

The  boundary  conditions  for  the  problems  are  met  by  the 
routine.   The  left  boundary  condition  assumes  that  the  slope 
is  zero  across  the  streamer  centerline.   So  here,  the  new 
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value  at  n(2,2)  is  assigned  to  n(l,2)  before  the  switch. 
The  right  boundary  assumes  that  n  is  zero  at  the  wall,  that 
is  station  101.   This  is  met  by  limiting  the  maximum  radius 
value  to  be  one  less  than  "RADIUS"  (called  "TEMRAD") .   The 
routine  procedes  to  "TEMRAD,"  then  the  last  point  seen  by 
the  calculation  algorithm  is  zero.   Similarly ,. the  tempera- 
ture of  the  right  boundary  is  set  to  1.0. 

To  compute  the  conductance  parameter,  gamma,  that  would 
result  from  any  particular  n  distribution,  call  "CNDUCT" 
after  assigning  the  desired  n  distribution  to  some  "VALX." 
The  single  value  "NN"  will  be  returned,  which  is  gamma. 

It  is  hoped  that  the  program  is  sufficiently  self- 
explanatory.   All  variable,  constant  and  subroutine  names 
are  defined  and  documentation  comments  are  provided. 
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APPENDIX  B 
HP-41  PROGRAM 

In  the  course  of  developing  the  results  of  Chapter  III, 
it  is  necessary  to  perform  a  large  number  of  simple  but 
tedious  calculations.   The  following  HP-41CV  program  proved 
helpful.   It  should  be  easily  adaptable  to  any  small  pro- 
grammable calculator. 

2  3 

The  program  prompts  the  user  for  E/N  (V-cm  ) ,  v . /N  (cm  /s) , 

-3         3 
N  (cm   )  ,  a~  (m  /s)  and  v   (cm/s) .   These  data  are  available 

from  Figure  2.2. 

-3         3 

The  program  then  displays  E  (V/m)  ,  N   (m   )  ,  a.      (m  /s )  , 

DTA,  BETAL  (laminar),  BETAT  (turbulent),  BETAAL  (laminar), 
BETAAT  (turbulent)  and  DT  in  true  time  in  microsecs .   The 
minimum  BETA  values  computed  are  consistent  with  the 
stability  criterion,  Equation  (2.17).   It  is  good  practice 
to  add  about  1%  to  these  values  to  ensure  a  smooth  running 
calculation . 
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01 

LBL  [Name  Program] 

37 

STO0  8 

73 

3.199E- 

02 

E/N? 

38 

RCL0  4 

74 

* 

03 

PROMPT 

39 

ENTER 

75 

ST013 

04 

STO01 

40 

RCL0  7 

76 

RCL10 

05 

V/N? 

41 

* 

77 

ENTER 

06 

PROMPT 

42 

i/x 

78 

RCL12 

07 

STO0  2 
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STO09 

79 

/ 

08 

N? 

44 

RCL09 

80 

ST014 

09 

PROMPT 

45 
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81 

200 

10 

STO0  3 

46 

10 

82 

ENTER 

11 

ALPHA2? 

47 
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83 

RCL14 

12 

PROMPT 

48 

STO10 

84 

* 

13 

STO0  4 

49 

1.6E-19 

85 

ST015 

14 

Vd? 

50 

ENTER 

86 

RCL0  8 

15 

PROMPT 

51 

RCL0  7 

87 

ENTER 

16 

STO0  5 

52 

* 

88 

RCL15 

17 

RCL01 

53 

RCL0  5 

89 

* 

18 

ENTER 

54 

* 

90 

RCL13 

19 

RCL0  3 

55 

RCL0  6 

91 

/ 

20 

* 

56 

/ 

92 

ST016 

21 

100 

57 

100 

93 

RCL0  6 

22 

* 

58 

/ 

94 

E  = 

23 

STO0  6 

59 

STOll 

95 

ARCL  X 

24 

RCL0  2 

60 

2.223E5 

96 

AVIEW 

25 

ENTER 

61 

ENTER 

97 

STOP 

26 

RCL0  3 

62 

RCL11 

98 

RCL0  7 

27 
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63 

/ 

99 

NO  = 

28 

RCL0  4 
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29 

/ 

65 
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AVIEW 
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STOP 
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RCL0  4 
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ALPHA2 
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RCL0  6 

106 

AVIEW 

35 

8E-3 

71 

X/*2 

107 

STOP 

36 

* 

72 

* 

108 

RCL14 
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109 

DTA  = 

110 

ARCL  X 

111 

AVIEW 

112 

STOP 

113 

RCL0  8 

114 

BETAL  = 

115 

ARCL  X 

116 

AVIEW 

117 

STOP 

118 

RCL16 

119 

BETAT  = 

120 

ARCL  X 

121 

AVIEW 

122 

STOP 

123 

RCL13 

124 

BETAAL  = 

125 

ARCL  X 

126 

AVIEW 

127 

STOP 

128 

RCL15 

129 

BETAAT  = 

130 

ARCL  X 

131 

AVIEW 
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STOP 

133 

RCL10 

134 

ENTER 

135 
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DTT  = 

138 

ARCL  X 
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AVIEW 

140 

END 
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